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ABSTRACT 

Young stellar systems orbiting in the potential of their birth cluster can accrete from the dense molecular 
interstellar medium during the period between the star's birth and the dispersal of the cluster's gas. Over 
this time, which may span several Myr, the amount of material accreted can rival the amount in the initial 
protoplanetary disk; the potential importance of this 'tail-end' accretion for planet formation was recently 
highlighted by Throop & Bally (2008). While accretion onto a point mass is successfully modeled by the 
classical Bondi-Hoyle-Lyttleton solutions, the more complicated case of accretion onto a star-disk system defies 
analytic solution. In this paper we investigate via direct hydrodynamic simulations the accretion of dense 
interstellar material onto a star with an associated gaseous protoplanetary disk. We discuss the changes to 
the structure of the accretion flow caused by the disk, and vice versa. We find that immersion in a dense 
accretion flow can redistribute disk material such that outer disk migrates inwards, increasing the inner disk 
surface density and reducing the outer radius. The accretion flow also triggers the development of spiral density 
features, and changes to the disk inclination. The mean accretion rate onto the star remains roughly the same 
with and without the presence of a disk. We discuss the potential impact of this process on planet formation, 
including the possibility of triggered gravitational instability; inclination differences between the disk and the 
star; and the appearance of spiral structure in a gravitationally stable system. 

Subject headings: accretion, accretion disks — planetary systems: formation — planetary systems: protoplan- 
etary disks 



1. INTRODUCTION 

Generally, stars form in clusters (e.g. Lada & Lada 2003), 
and our current understanding of star formation is beginning 
to incorporate the effects of the cluster environment on the 
process rather than treating star formation as a one-system 
event occurring in isolation. A complete picture of planet for- 
mation may require a similar awareness of the surroundings; 
the possible influence of the cluster on planet formation is a 
logical consequence of clustered star formation. Protoplane- 
tary disks have lifetimes of ~ 1-10 Myr without significant 
external UV irradiation (e.g. Currie et al. 2007, 2009). In the 
Orion Nebula Cluster, 80% of stars have disks (Smith et al. 
2005), and even in environments like the Trapezium clus- 
ter over 10% of stars may have disks containing a minimum 
mass solar nebula (MMSN, ~ 0.01 M Q ) of material (Mann 
& Williams 2009). These disk lifetimes are comparable to or 
longer than the crossing time of a typical star-forming region, 
and thus the cluster environment may affect planet formation 
or young planetary systems. 

Disruptive possibilities include the truncation of disks by 
interactions with other cluster members (Clarke & Pringle 
1993; Heller 1995; Boffin et al. 1998) or dynamically alter- 
ing an existing planetary system (Laughlin & Adams 1998; 
Zakamska & Tremaine 2004; Malmberg & Davies 2009) \ 
and photo-evaporation which can remove disks on 10 5 -10 7 
yr timescales (Hollenbach et al. 1994; Johnstone et al. 1998). 

Electronic address: nickolasl@gmail.com 

1 Although the rates for these interactions are probably low compared to 
formation timescales (Adams et al. 2006). 



More constructive effects may include triggering planetesi- 
mal formation by external UV illumination (Throop & Bally 
2005), and adding material to the disk by accretion from the 
cluster's remnant ISM reservoir (Throop & Bally 2008, here- 
after TB08) 

The effect that we explore in this work is the accretion of 
ISM material by a star-disk system. Bondi-Hoyle-Lyttleton 
accretion by a star moving through a background medium 
may play a fundamental role in determining the final stellar 
mass (Bonnell et al. 1997, 2001), an idea that has seen a large 
amount of study over the last decade. A less-explored ques- 
tion is the importance of accretion after a star has already as- 
sembled the majority of its mass, and the cluster is dominated 
by the gravitational dynamics of the stars rather than the po- 
tential and bulk velocities of the molecular gas. 

TB08 performed simulations of this 'tail-end' accretion, 
finding that as stars orbit through the depleting molecular gas 
in a cluster they can accrete around 1 MMSN per 1 Myr across 
a range of standard star-forming regions. The observational 
constraints on disk lifetimes allow for this accretion to take 
place while a primordial protoplanetary disk is still present. 
A plausible example of the interaction between a disk and the 
ISM is found around the star HD 61005 (Hines et al. 2007), 
displaying what appears to be a wind-swept disk. The TB08 
study highlighted the potential importance of this process and 
suggested several possible measurable effects, including al- 
tering the structure of the protoplanetary disk, changing its 
orientation, or changing its composition. 

It is these ideas that we explore here via direct hydrody- 
namic simulations. TB08 treated accretion by an analytic pre- 
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scription applied post-simulation to n-body experiments, and 
here we present the results of exploratory simulations involv- 
ing a protoplanetary disk that encounters a dense ISM at su- 
personic speeds. In TB08 the accretion is treated as occurring 
from a mean background density over ~ 10 6 yr; the varia- 
tion we consider is a star-disk system encountering a dense 
filament or clump, with a higher accretion rate over a shorter 
amount of time. After describing our setup and presenting the 
results of our experiments, we discuss their implications for 
theories of planet formation. 

2. PHYSICAL AND NUMERICAL MODEL 

2.1. Bondi-Hoyle-Lyttleton Accretion 

Accretion onto a point mass moving with some velocity 
through an ambient medium is a well-studied question, and 
the classical solutions developed by Bondi, Hoyle and Lyttle- 
ton (BHL) (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; 
Bondi 1952) are widely used as a starting point (and often as 
an end point) when considering accretion problems. Many 
numerical tests have added to our understanding, offering 
verification and modification to the analytic accretion rates 
(Shima et al. 1985; Ruffert 1994b; Ruffert & Arnett 1994; 
Ruffert 1994a, 1995, 1996, 1999) and considering compli- 
cations like the effects of radiation (Edgar & Clarke 2004), 
vorticity (Krumholz et al. 2005), and turbulence (Krumholz 
et al. 2006). A good review of the historical development and 
numerical testing of BHL accretion theory is found in Edgar 
(2004). 

The general picture of BHL accretion is of a point mass in 
motion through a uniform, infinite medium. Ambient mate- 
rial is focused into a wake or accretion column behind the 
accretor, where colliding streamlines cancel out velocities 
that aren't parallel to the body's motion. Material passing 
within some critical radius becomes bound and is energeti- 
cally doomed to flow up the wake and be accreted by the star. 
Given the density and sound speed of the ambient medium, 
and the velocity and mass of the accretor, the critical radius 
can be calculated and thus the accretion rate. The analytic ex- 
pression for the accretion rate varies by a few factors of 2 in 
the literature, and we briefly explain our choices here. 

The critical impact parameter of a test particle moving at 
velocity v past a body of mass M can be found by simply 
equating the potential and kinetic energies, and one finds the 
Hoyle-Lyttleton radius R HL = 2GM/v 2 (Hoyle & Lyttleton 
1939). For highly supersonic motion, accretion is firmly in the 
Hoyle-Lyttleton limit and this would be a fine approximation 
to the accretion radius for this work. For ease of comparison 
to recent numerical results, we will use as our accretion ra- 
dius the interpolation between the supersonic and stationary 
accretion regimes proposed by Bondi (1952) with the factor 
of 2 modification suggested by Shima et al. (1985), 



and an associated accretion rate 

M BH = irR 2 A p(\ 2 c 2 s +v 2 )V 2 . (2) 

The constant A is of order unity and depends on the equation 
of state of the gas; for the isothermal limit, A = e 3/2 /4 ~ 1 . 120 
(Bondi 1952). The mass density and sound speed of the am- 
bient medium are p and c s , respectively. 

2.2. Parameters and Initial Conditions 



Throughout this work, we take the mass of the accretor to 
be M = 1 M Q , and its velocity through the medium v = 3 km 
s" 1 . We assume the medium has a mean molecular weight 
/j, = 2.3 and a sound speed c s = 0.3 km s" 1 , thus a temperature 
~ 25 K and a Mach number of 10. The number density of 
the medium is set at n = 5 x 10 6 cm" 3 , yielding a mass density 
/9=1.92xl0~ 17 g cm" 3 . Plugging these values into equations 
1 and 2 gives us our reference accretion radius and rate and 
the natural timescale for the problem, 

Ra = 195 AU 
M BH = 2.47 x 10" 6 M yr -1 
t bh = R a /c s ^ 3000 yr. (3) 

These ISM density values are quite high, and are more asso- 
ciated with the regions surrounding dense cores (which can 
reach mean hydrogen densities in excess of 10 7 cm" 3 at 0.05 
pc from the core, e.g. Mezger et al. 1990) or over-dense ridges 
(e.g. a 0.5 pc ridge with n > 10 6 cm" 3 near the BN/KL region, 
Kawamura et al. 2002), rather than the ambient background. 
Thus our simulations are only directly applicable in the dens- 
est regions of molecular clouds. This ISM density is picked 
for computational reasons discussed below, rather than purely 
physical motivation. 

Our simulations model a gaseous protoplanetary disk im- 
mersed in the ISM flow. The disk initially contains 0.01 M Q , 
extending from 10 - 100 AU with a surface density profile 
S oc r" 3 ' 2 . For reference, this is comparable to the density of 
a classical MMSN (Weidenschilling 1977), and would give a 
surface density ~ 1000 gm cm" 2 at 1 AU were the disk to ex- 
tend in that far. The disk and the inflowing material are both 
isothermal at the same temperature. 

We carried out four simulations: accretion from a disk with 
no ambient medium; accretion onto a point mass (BHL accre- 
tion); BHL accretion with a disk perpendicular to the inflow; 
and BHL accretion with the disk inclined 30° to the inflow. 
The first two experiments serve to set baseline accretion rates 
onto the star from the ISM and the disk in isolation, while 
the experiments with both inflow and a disk explore how the 
disk's presence alters the ISM accretion flow and vice versa. 
We run each simulation until t « 5000 yr, during which the 
star or star-disk is immersed in the accretion flow for 1 t B h 
and evolves for a further 0.5 t B h after exiting the ISM. This 
is sufficient time for a quasi-stable accretion state to be set 
up. We will refer to time in years as t , and time in units of t B h 
from as r. The zero point for each timescale is set to when the 
ISM flow first reaches the star. Each system evolves in effec- 
tive isolation for t ~ 525 yr leading up to t = 0, which allows 
the disks to equilibrate and gives the ISM time to traverse 
from the inflow boundary to the star's position. The inflow 
boundary is turned off such that after being immersed in the 
flow for 1 t B h the star exits the dense medium, the remain- 
ing gas is accreted or leaves the domain, and the disk settles 
down. We can then see the effect on a protoplanetary disk of 
a passage through a dense clump or filament in a molecular 
cloud. 

In TB08, the motivating work for this study, the stars or- 
bited and accreted over a period of 4 Myr, accreting from a 
Plummer sphere gas distribution with a peak density ranging 
from 10 4 — 10 6 cm" 3 depending on the total cluster mass. Our 
background density is five times higher than the nominal peak 
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density 2 in an Orion molecular cloud analog, 10 6 cm -3 , thus 
our picture of a collision with a prestellar core or filament. 
This choice is motivated by computational necessity rather 
than just physical arguments, as such collisions will be very 
rare in all but the densest protoclusters. A simple estimate 3 of 
the encounter rate between a protostar and a collapsing core 
yields encounter times of at least ~ 0.7 Myr. While this is 
short enough to allow for a handful of such encounters over 
the cluster lifetime assuming continuing star formation for a 
few Myr, it is longer than the likely collapse time for dense 
prestellar cores (a few xl0 4 -10 5 yr, e.g. Ward-Thompson 
et al. 2007), which is the more relevant time constraint. As 
described more fully in section 2.4, we have chosen to mini- 
mize the density disparity between the disk and ISM, pushing 
the disk mass to low values and the ISM density to unrealis- 
tically high values. Additionally, we simulate the interaction 
only long enough to set up a quasi-stable accretion flow be- 
fore ending the simulations, so that the disk is left as unaf- 
fected as possible by the numerical viscosity inherent in any 
computational scheme. These are the constraints that moti- 
vate our scenario of an encounter with a dense clump or core 
rather than the background ISM as in TB08. The time-scale 
arguments above suggest that these specific simulation param- 
eters are directly applicable to only a few percent of stars in 
an Orion-like environment. However, as we argue in the next 
section, the main effect of the ISM accretion on the disk is 
not impulsive but rather a function of the total mass accreted, 
and so the results of these over-dense simulations should be 
broadly applicable to the original scenario in TB08. 

2.3. Expected Effects 

There are several effects that one can imagine occurring 
when a gas disk is subjected to what is effectively a super- 
sonic wind, including modifications to the disk structure and 
the morphology of the accretion column. Here we discuss a 
few of them and estimate how likely they are to play a role 
given our chosen parameters. 

Erosion by ram pressure stripping: Depending on the pa- 
rameters of the ISM, disk density, stellar mass and velocity, 
the wind may truncate the disk at some radius. Ram pressure 
stripping is frequently considered in studies of galaxies in- 
teracting with the intra cluster medium, a situation for which 
several analytic estimates of the truncation radius have been 
developed (Gunn & Gott 1972; Mori & Burkert 2000). More 
relevant to this work, the question of supernovae shells in- 
teracting with protostellar disks was discussed by Chevalier 
(2000). As in these works, we compare the ram pressure of 
the ISM flow to the gravitational restoring force keeping the 
disk associated with the star. The radius at which these are 
comparable should give some indication of the importance of 
ram pressure stripping to our 100 AU disk, which we assume 
here to be face-on to the ISM flow. 

The gravitational force per unit area (the gravitational 'pres- 
sure' to which we can compare the ISM ram pressure) on 

2 A Plummer sphere with 3000 Mq of gas matching Orion's current core 
radius of ~ 0.2 pc (Hillenbrand & Hartmann 1998) has a central density 
~ 10 6 cm 4 . 

3 For example, taking the m = 1 Mq protostars or cores to have a colli- 
sion radius r = 0.01 pc, number density n+ = 10 4 pc~ 3 , and velocity disper- 
sion a = 2km s~' (appropriate for the Trapezium), the standard collision time 
(Binney & Tremaine 2008) is Q 1 ,, = 4 v / ¥«*cr(r 2 + Gmr/cr 2 ) = 6.6 X 10 5 yr. 
The scaling is linear with number density, so a cluster with = 10 3 would 
have an encounter rate in excess of 6 Myr. 



a ring of radius r is given by P grav = GMYi{f)jr 2 . For our 
disk with E(r) « 1000(r/AU)" 3 / 2 g cm" 2 , this is P grav « 
600(r/AU)~ 7//2 dyne cm" 2 . The ram pressure of the ISM flow 
is given by Pjsm = pv 2 w 1.7 x 10~ 6 dyne cm" 2 . The radius 
at which P ISM > P gm v is approximately 275 AU, larger than 
the initial disk radius of 100 AU. Bearing in mind the order- 
of-magnitude nature of this calculation, especially when es- 
timating Pgrav, it would not be surprising to see ram pres- 
sure affecting the extremities of the disk, but it seems clear 
that its effect should be negligible for the inner regions of the 
disk. Note that outside of 75 AU the gas pressure in the disk 
Pdisk = c 2 p < 1.7 x 10" 7 dyne cm -2 , one order of magnitude 
less than Pism- Flattening, but not unbinding, of the disk at 
the outer edge is to be expected. 

Momentum stripping: The ram pressure stripping calcula- 
tion assumes that a quasi-equilibrium exists between the pres- 
sure of the ISM wind and the gravity in the star-disk system, 
and thus that nothing changes on roughly the orbital timescale 
of the disk. A potentially shorter-timescale effect (shorter 
than, say, the the orbital period at the outer disk radius) is 
the addition of momentum perpendicular to the disk rotation. 
If the momentum addition raises the disk material to escape 
velocities, the disk will be eroded from the outer edge (Cheva- 
lier 2000). 

The escape velocity of disk material at radius r is v esc = 
(2GM I 'r) l l 2 . Since the material is moving at the Keplerian ve- 
locity already, assuming again a face-on disk we require a ve- 
locity gain perpendicular to the disk plane 5v z (r) = {GMj r) 1 / 2 
on a timescale shorter than the orbital period to unbind disk 
material at radius r. The maximal change in velocity would 
occur if the disk swept up much more than its own mass 
(which it won't with our parameters), in which case the disk 
material would move with the ISM velocity of 3 km s" 1 . This 
is less than the required velocity change even at 100 AU, so 
the disk is very robust to erosion by absorbing momentum 
from the ISM. 

Disk redistribution: While absorbing momentum from the 
ISM will not unbind our disk at any radius, it may have an 
effect on the radial distribution of material. The disk mass as- 
sociated with the annulus of infinitesimal radius dr at radius 
r is rridisk = 2irrT,(r)dr, while the swept-up ISM material af- 
ter some time t is m iSM = 2-irrpvtdr. When the ISM material 
is absorbed, the specific angular momentum in that annulus 
drops by a factor of S(r)/[S(r) + pvf]. Assuming that the an- 
nulus adjusts to the radius whose Keplerian angular momen- 
tum reflects this new value, the new radius f can be found 
from r l l 2 w r^ 2 T,(r)/[T,(r) + pvt]. 

This rough calculation ignores any interaction with neigh- 
boring annuli, and assumes a smooth addition of material to 
the disk. Those caveats in mind, after accreting for a time 
t = tbh, material at 50 AU will have migrated to 48 AU, while 
the outer edge of the disk will have moved from 100 AU to 
90 AU. The trend of this effect is for the disk to become more 
compact, with the outer edges migrating proportionally fur- 
ther in. 

Summarizing the likely effects of the ISM flow on the disk: 
ram pressure will only noticeably alter the outskirts of the 
disk. Mass and momentum loading will not unbind the disk in 
any way, but will redistribute material inwards, and this effect 
will increase with radius. We see that even at the very high 
ISM densities we are simulating, the ram pressure is not the 
dominant effect. Rather, the total mass that is absorbed (over 
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timescales longer than the orbital period of the disk) domi- 
nates the disk evolution. At more realistic ISM densities the 
ram pressure effects, already less important than mass load- 
ing, will be further diminished. The total mass accreted by 
our systems over their 3000 yr immersion in the dense ISM, 
roughly M^z/f ~7.5 x 10~ 3 M©, is less than that accreted by 
the average solar-mass star in TB08, and thus the disk redis- 
tribution should be comparable to that seen here. 

Alteration of the accretion flow: We estimated above the 
unimportance of the ISM ram pressure in influencing the disk 
evolution; conversely, this suggests that the disk will alter the 
accretion flow compared to the diskless case. Since the disk 
radius is r^k ~ Ra/2 and is centered on the axis of the ac- 
cretion flow, one-quarter of the accretion flow is intercepted 
by the disk. The intercepted material is on streamlines that 
focus to the accretion axis close in to the star. Without these 
converging streamlines compressing the wake, the accretion 
structure close to the star should be less collimated in the 
presence of a disk, as the final inflow occurs in the 'lee' of 
the disk. 

2.4. The Code 

We used the smoothed particle hydrodynamics (SPH) code 
GADGET-2 (Springel 2005) to perform our 3D hydrody- 
namic simulations. The code has been modified to suit this 
problem in the following ways: 

Sink particles: The accretor is treated as a sink particle 
(Bate et al. 1995), accreting all particles that pass within its 
accretion radius and absorbing their linear momentum. By 
removing particles that would otherwise have a tiny timestep 
we are able to follow the evolution of the system for a longer 
time. We apply no boundary conditions at the accretion ra- 
dius, which we set at R acc = 9.75 AU = 0.05 Ra\ this value is a 
compromise between computational efficiency and the reality 
of a nearly-point mass accretor. The comprehensive series of 
BHL accretion simulations by Ruffert showed that with ac- 
cretors larger than this, the character of the flow is fundamen- 
tally changed as accretion is onto the leading edge of the body 
rather than along the accretion column. 

Artificial viscosity: SPH codes rely on an artificial viscos- 
ity to capture shocks and prevent particle interpenetration. 
In place of the GADGET-2 default viscosity, we have im- 
plemented the viscosity suggested by Morris & Monaghan 
(1997). Here each particle has its own viscous parameter a, 
which evolves according to 

d = -(a-0.1)/T a + max(-|V-v|,0)(2-a). 

The first term decays the viscosity down to the minimum 
value 0.1 over the timescale r a , while the second term is 
meant to detect shocks and cause the viscosity to grow. Each 
particle's individual a then figures in the standard viscosity 
formulation, rather than a fixed value for all particles. This 
reduction in the numerical viscosity away from shocks gives 
a less dissipative treatment of smooth flows and quiescently 
rotating disks. 

Inflowing and outflowing boundary conditions: A critical 
element of our study is the inflow of ambient ISM material 
into the simulation domain. We have implemented bound- 
ary conditions such that the star sits in a cylinder with height 
and diameter 5 Ra- The star is placed on the cylinder's axis, 
2.5 Ra from both the inflow and outflow boundaries. The 
ISM flow is introduced at the inflow boundary and consists 
of randomly placed SPH particles at the chosen velocity and 
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Figure 1. Top: column density plot of Mach 10 BHL accretion, showing 
only material 2.5 times denser than the background density. The ISM flows 
in from the right and exits left. The accretor is shown as a gray circle. Bottom: 
black arrows indicate the velocity field of the flow, showing direction only. 
Solid grayscale streamlines are the analytic flow field of the Hoyle-Lyttleton 
solution. Dark gray streamlines are interior to the critical impact parameter, 
while unbound streamlines are in light gray. The stagnation point of the sim- 
ulated accretion column matches the analytic value, though the streamlines 
differ somewhat near the accretion axis. 

mass. A boundary layer of particles that experience no ac- 
celeration shepherds the ISM flow at the outer edge of the 
cylinder. SPH particles that leave the domain through the out- 
flow edge are removed from the calculation. These are not 
periodic boundary conditions; the ISM material introduced at 
the inflow boundary is 'fresh' . 

This calculation is challenging in that there is a large mis- 
match between the densities of the ISM and the disk, which 
we have ameliorated as much as possible by choosing a dense 
ISM and a low-mass disk. SPH codes use a Lagrangian 
scheme in which the fluid is represented by discrete tracer 
particles, which are interpolated over a density-dependent 
smoothing length to obtain continuum values throughout the 
domain. Choosing equal mass ISM and disk particles yields 
an undesirably large mismatch in smoothing lengths, and set- 
ting the masses such that the smoothing lengths are equal 
gives unacceptably disparate masses. By setting the disk par- 
ticle mass to be 8 times the ISM particle mass we avoid the 
problems associated with large particle mass ratios, while still 
having the ISM particle smoothing lengths set such that the 
ISM flow comfortably resolves the protoplanetary disk. The 
disk is initially modeled with ~ 2.5 x 10 particles; the total 
particle number is constantly varying, but when the cylinder 
is full there are ~ 5 x 10 6 ISM particles. 

3. NUMERICAL RESULTS 

3.1. BHL Accretion With No Disk 

We begin by examining the morphology of the accretion 
flow with no disk present. Our isothermal inflow is Mach 10 
and our accretor is 0.05 Ra in radius, a setup which lies be- 
tween models HS and HM in Ruffert (1996), to which we 
can compare our results. Figure 1 shows the column density 4 
and velocity field for this run once the accretion flow is estab- 
lished. In this, and all surface density plots in the paper, we 
plot only the contribution to the column density from particles 
with density greater than 2.5 times the background. This ef- 

4 We used the progTam SPLASH (Price 2007) to project the SPH particle 
positions into column density values and velocity fields for this and similar 
plots. 
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Figure 2. Accretion rates onto the star from the two sources in the simulations; black lines are accretion from the ISM, while gray lines are accretion from the 
disk. Arrows at the left of each panel show the mean ISM accretion rate for times 0.5 < t/t bh < 1.0. The time scale at the top of the panels is solid when the 
ISM flow is on, and dashed when the star has left the flow; the star enters the dense ISM at at r = 0, and remains in it until r = tbh, and we ran the simulations 
for an additional 0.5tbh after the end of the ISM flow has passed by the star. Top left: accretion from the ISM with no disk. Bottom left: baseline accretion 
rate from the disk with no ISM. Top right: both accretion rates for the case with the disk face-on to the inflow. Bottom right: both accretion rates with the disk 
inclined 30° to the inflow. 



fectively isolates the shocked accretion column from the sur- 
rounding ISM, removing the arbitrary scaling associated with 
the size of the accretion domain as well as removing some ar- 
tifacts from the boundary layer of particles. The simulation 
domain is larger than the region shown; we show here just the 
accretion column within 2 Ra of the star, while the simulation 
extends 2.5 Ra in radius, upstream, and downstream from the 
star. 

The column density plot appears largely as one would ex- 
pect. With high Mach number flows the accretion column has 
a small opening angle, tightly compressed behind the accretor 
and widening as the distance from the accretor increases. We 
do see larger breaks from axisymmetry in our run than Ruffert 
did, an effect that also manifests itself in our accretion rates. 
A possible cause of this is the way we introduce the accretion 
flow, randomly seeded in space rather than on a fixed grid. 
The random placement of our SPH particles, while maintain- 
ing an overall mass inflow rate that matches our chosen value, 
allows for random pockets of over- and under-density that 
do not smooth themselves out before reaching the accretion 
shock. These may act to seed instability in the accretion col- 
umn. 

The velocity field shows the fluid direction only, not the 
magnitude, overlaid on the analytic streamlines of the Hoyle- 
Lyttleton solution from Bisnovatyi-Kogan et al. (1979). The 
darker streamlines are expected to remain bound to the ac- 
cretor, while lighter streamlines should escape. The position 
of the stagnation point in our simulation is variable in time, 
but always remains in the vicinity of the expected location. 
While the character of the ballistic streamlines matches the 
simulation, the detailed shape is different, especially as the 
flow approaches the accretion axis. 



In figure 2 we show the accretion rates for all four simu- 
lations. The upper left panel shows the BHL accretion case. 
There is no accretion until the inflow turns on and reaches 
the star at t = 0. The accretion rate ramps up to an oscil- 
lating, quasi-steady state at t ~ 0.5 t bh . The large spike in 
the accretion rate at r ~ 0.3 tbh is from when a clump that 
formed as the edge of the accretion flow converged down- 
stream, stagnated, and finally accreted back onto the star. We 
believe this to be an artifact of the way the star effectively 
enters the dense region abruptly from a zero-density region, 
rather than beginning immersed in the ambient medium. The 
increase in the accretion rate at r <~ 1.25 tbh is a result of 
the star leaving the dense ISM. Gas velocities near the stag- 
nation point are subsonic, and near the end of the accretion 
period a large blob of gas forms and drifts slowly away from 
the star and out of the domain. When the dense ISM ends and 
there is no pressure from the upstream direction, this gas flows 
back toward the star and is accreting as the simulation ends. 
We measure the mean accretion rate onto the star between 
0.5 < t/tbh < 1, i-e. between when the accretion column is 
fully established and before the stars leaves the dense ISM, as 
< M >= O.IIMbh- This agrees within ~ 6% of the roughly 
equivalent case in Ruffert (1996) of < M >= 0.82M BH , though 
we have more variation about the mean than in his simula- 
tions. 

3.2. Disk Evolution With No Inflow 

In order to get a baseline accretion rate from the disk with 
no ISM accretion, we ran the disk in isolation for the same 
amount of time. The accretion rate is shown in the lower left 
panel of figure 2. The accretion rate is marked by a high initial 
value as the inner edge of the disk drains onto the star, before 
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Figure 3. Top: azimuthally-averaged surface density profile for the disk 
evolving in isolation. Middle: evolution of the disk face-on to the inflow. The 
gray levels are as in the top panel. Solid lines are the total surface density, 
and dashed lines are the contribution to the disk from the ISM. Bottom: lines 
are as for the middle panel, but for the case with the disk inclined 30° to the 
inflow. 

settling down to a low value of ~ 1 .5 x 10~ 7 M Q yr" 1 . 

The top panel of figure 3 shows the evolution of the 
azimuthally-averaged surface density profile for the isolated 
disk case. The lightest line shows the profile just after the sim- 
ulation began, and further lines are at values of t bh / 2 starting 
at t = (recall that this is the time when the accretion flow en- 
counters the disk). The wiggle in the density profile at t = 
is the final transient of the disk equilibration. Throughout the 
run the inner edge of the disk drains onto the accretor. This 
is due predominantly to the numerical dissipation from the ar- 
tificial viscosity, with a small contribution from the lack of 
gas pressure support from SPH particles starward of the in- 
ner edge. Boundary conditions at the accretor might prevent 
this drainage, but rather than risk introducing other effects we 
allow the disk to evolve under the influence of artificial vis- 
cosity. There is some spreading at the outer edge as well, but 
the density profile over the bulk of the disk remains almost 
unchanged. Throughout the simulation the disk is stable to 
gravitational instabilities, and no spiral features develop. 

3.3. BHL Accretion With Disk Inclined 0° To Inflow 

We turn now to the runs including both a disk and inflow, 
and consider the case of a disk with angular momentum par- 
allel to the motion of the star through the ISM. With this ar- 
rangement the disk is maximally absorbing the column of ra- 
dius Ra that the star is drilling through the ISM. Recall that 
our choice of parameters has Rdisk/RA ~ 0.5. This means that 
material with impact parameters < 0.5Ra will be intercepted 
by the disk rather than proceeding to collide downstream and 
accrete up the wake. The linear momentum of this material is 
absorbed by the disk, leading to changes in the disk density 
profile and stripping of the outer regions. 

Figure 4 shows the morphology of the disk and the accre- 
tion column as it develops. We plot the projected surface den- 
sity through the simulation domain on a logarithmic scale. 
The top row shows the disk just before it begins interacting 
with the leading edge of the ISM flow. In the second row, at 



time r = 0.25tbh, the accretion flow is beginning to assemble 
downstream of the star, and the face-on projection shows that 
the disk structure is changing under the influence of the accre- 
tion flow. By r = 0.5t bh the accretion column is established, 
and much like the pure BHL accretion case accretion proceeds 
up the wake and onto the accretor. The main morphological 
difference is the looser collimation of the accretion column 
close to the star. With the disk intercepting the inner part of 
the upstream flow, the wake is accreting in the lee of the disk 
and the strong shock from material with a small impact pa- 
rameter is unavailable to constrain the flow. This results in 
the wispy appearance of the accretion stream close to the star. 
The stagnation point of the downstream flow, as in the diskless 
case, is in the vicinity of 1 R^ from the star. A slight drift in 
the position of the star in the direction of the ISM flow is also 
seen. In our physical picture, this is due to the star slowing 
down as it interacts gravitationally with the wake of material 
behind it. 

The upper right panel of figure 2 shows the accretion rates 
onto the star from the disk material and the ISM. The growth 
of the ISM accretion rate has a different shape in the presence 
of a disk compared to the BHL case. In the diskless case, the 
early accretion as the wake is being established is due to low 
impact parameter material with colliding streamlines close to 
the star. With a disk present, this material is intercepted and 
the final quasi-steady accretion flow takes longer to develop. 
The low accretion rate in this case from < t/t bh < 0.25 is 
from ISM material with impact parameter less than the accre- 
tion radius of our sink particle, which accretes from the up- 
stream side of the star. After r ~ 0.25tbh the accretion flow 
generated by larger impact parameter material that avoided 
the disk begins to establish itself, and the ramp-up to a quasi- 
steady state begins. From r = 0.5t B h until the ISM flow stops 
the accretion rate oscillates unstably around a mean value of 
2.25 x 10" 6 M yr" 1 (just under the reference accretion rate 
of equation 3), though there does appear to be an increasing 
trend over this interval. 

Accretion from the disk is largely the same as the baseline 
rate, with some notable spikes that coincide with peaks in the 
ISM accretion rate. These episodes of increased accretion oc- 
cur when a dense stream of ISM material, unconstrained ra- 
dially as it approaches the star, impacts the disk at a radius of 
~ 20 AU. The addition of this low angular momentum mate- 
rial to the disk causes an accretion event. Note that the tim- 
ing and details of this process are dependent on the numerical 
accretion radius. It is possible that this material would have 
adjusted to a smaller radius inside of our 10 AU accretion 
radius, and likewise this same effect (ISM accretion streams 
striking the disk) could occur at smaller radii and scales than 
we follow in these simulations. 

After the flow is turned off, the remnants of the accretion 
stream either escape from the system or accrete onto the star, 
leaving the disk to equilibrate. In the accretion rate plotted in 
figure 2, the accretion of the final lingering ISM wisps is seen 
as a non-zero ISM accretion rate, but this is not affecting the 
remaining disk in any significant way. In the final panel of fig- 
ure 4, this low-density ISM material is visible in the edge-on 
view, with the denser material on the left edge moving away 
from the star toward the outflow boundary. The face-on view 
shows no remaining spiral structure, as the disk has settled 
down to its final state. 

The middle panel of figure 3 shows the surface density of 
the disk during this run. Because the disk is changing shape 



Figure 4. Snapshots from the simulation with the disk face on to the inflow, showing the column density on a logarithmic scale. The left column is a projection 
orthogonal to the ISM flow vector, which is from right to left. The right column is the projection along the direction of the ISM flow. Time increases downward, 
as labeled. Only material denser than 2.5 times the background density contributes to the column density. An animation of this simulation can be found at 

http: //origins. Colorado. edu/~moeckel /BHLpaper /. 



and perhaps orientation due to the accretion flow, we define 
the disk for this plot as follows: considering only material 
within 120 AU of the star, the angular momentum of all re- 
maining disk particles are summed. Then all ISM particles 
with angular momentum orientations within 7r/8 of this vec- 
tor are found, and included in the surface density calculation. 
This procedure is effective in removing the accretion stream 
and isolating what one would visually deem part of the disk. 
The solid lines in the plot show the total disk surface density, 
i.e. both the ISM and the original disk particles. The dashed 
lines show only the contribution to the disk from the absorbed 



ISM. 

The two most obvious features of the disk evolution are the 
inward migration of the outer edge, and the upward trend in 
the surface density profile. While it is initially tempting to 
attribute these effects to truncation and addition of swept up 
material, respectively, the dominant cause of both is the redis- 
tribution of disk material, with a contribution to the increased 
surface density from adding ISM material to the disk. To il- 
lustrate this, compare the total surface density profile at the 
initial time to the profile at t = 0.5tbh in figure 3. At 50 
AU, the total disk surface density has increased from ^3 g 
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cm" 2 to ^7 g cm -2 , while the contribution from the ISM is 
just under 0.9 g cm" 2 ; the remaining increase is from mate- 
rial that was originally part of the disk, but moved inwards 
from larger radii. From t=1.0-1.5t bh the ISM contribution 
appears to decrease, but this just matches the viscous draining 
and spreading of the disk, as the inflow is no longer accreting 
onto it. The argument against truncation being the cause of 
the smaller outer radius is bolstered by the amount of original 
disk mass remaining. If the disk were truncated at ~ 60 AU, 
as might be guessed when looking at the plot, the remaining 
disk mass should be 67% of the initial amount. Instead, over 
90% remains, and the missing 10% is due to accretion rather 
than loss to the ISM. As estimated in section 2.3, ram pressure 
is ineffectual at actually stripping the disk, and inward migra- 
tion of the disk due to mass loading is occurring, affecting the 
outer radii more than the inner disk. The final disk has a shal- 
lower and denser profile than it started with, and depending 
on the radius ~ 10-25% of the disk material came from the 
ISM. 

The full situation is of course more complicated than our 
simple estimates suggested. While ram pressure cannot re- 
move the outer radii of the disk, it does have an early influ- 
ence on the disk morphology, as is clearly seen when viewing 
a movie of the simulation snapshots (and perhaps in the real 
universe, as in HD 61005; Hines et al. 2007). The accretion 
flow has a velocity component that is directed radially inward 
by the time it hits the disk. The outermost 10 AU or so to 
remain bound to the star-disk system, but they are deflected 
along the flow in toward the star, rejoining the disk plane in- 
ward of their original location and hastening the radial evolu- 
tion compared to our estimated rate. 

3.4. BHL Accretion With Disk Inclined 30° To Inflow 

At a bulk level, the accretion flow in the inclined case 
(shown in figure 5) is very similar to the face-on case, with 
a similarly shaped accretion column leading up to the star, 
and a stagnation point as expected near R A downstream of 
the star. The similarity is expected, since except for some 
minor geometrical effects from the disk inclination the situa- 
tions are nearly identical, with the inner half of the accretion 
radius intercepted by the disk. The accretion rates from the 
disk and the ISM, plotted in the lower right panel of figure 2, 
is also quite similar, with nearly identical average accretion 
rates from t=0.5-1tbh- Spikes in the ISM accretion rate are 
matched by spikes in the disk accretion rate, as the inner disk 
absorbs the low specific angular momentum material of the 
accretion column and is accreted. The same caveats about the 
size of our accretion radius apply to this case as the face-on 
simulation. 

While the accretion flow is largely the same between the 
two cases, one may expect the disk evolution to be different. 
The component of the Keplerian disk velocity in the direction 
of the inflow is greater than ~ 2 km s" 1 throughout the disk, 
and indeed greater than the inflow velocity of 3 km s" 1 inward 
of ~ 25 AU. While the absorption of the smooth inflow from 
the upstream direction should be averaged over the orbital 
period and negate this effect of the inclination, the effect of 
absorbing the dense and filamentary accretion column may 
depend on where it strikes the disk in a way not seen in the 
face-on case. Examining the surface density evolution in fig- 
ure 3 it appears that there may be some difference in the later 
evolution of the inner disk, but the difference is small enough 
that run-to-run variability could be the main cause. Future 



investigations will pursue this matter further. 

4. DISCUSSION AND CONCLUSIONS 

As discussed above, we were constrained by computational 
necessity to consider a disk interacting with a dense molec- 
ular ISM, corresponding perhaps to a filament or core in a 
star-forming region. The analysis of TB08 which motivated 
this work dealt with an ISM density reflecting the mean value 
of a region, much lower than this work (10 3 -10 6 cm" 3 ), and 
over longer timescales (Myr versus kyr). Because of the na- 
ture of the interactions, we believe that the main results of this 
work should carry over from the dense interactions we simu- 
lated to the lower density, longer timescale interactions from 
TB08. If ram pressure had been the dominant driver of the 
changes to the disk morphology, this would not be the case, 
so interactions with a much denser region or at much higher 
velocities will be qualitatively different than these. However, 
the changes were instead chiefly a result of the deposition 
of material with no azimuthal angular momentum onto the 
disk over timescales longer than (or comparable to) the or- 
bital timescale. These effects should remain broadly the same, 
whether the low specific angular momentum material is incor- 
porated into the disk over 10 or 1000 orbital periods. Like- 
wise, if the angular momentum is deposited episodically (as 
in TB08) versus in one interaction as we have modeled, the 
final effect should mostly be a function of the total mass de- 
posited. Thus the main result of this work, the redistribution 
of the outer disk material inward, should be robust over a wide 
range of ISM parameters. 

A potentially significant complication that we did not con- 
sider is that the ISM in a real star-forming cloud will not be 
smooth, but instead filamentary and clumpy. To some extent 
we have modeled this, as the star moved from a region of zero 
density into the dense region, and left after moving through 
it for ~ 3000 yr (equivalent to a core or filament ~ 0.01 pc 
across). However, the effect of substructure on smaller scales 
would be an interesting direction for future studies. Density 
or velocity structure on scales comparable to the accretion ra- 
dius could have an important effect on how the disk restruc- 
turing proceeds, introducing net angular momentum to the in- 
flow and perhaps changing the disk structure more dramat- 
ically than the almost totally-smooth ISM we have consid- 
ered. Simulations of density and velocity gradients in pure 
BHL accretion (Ruffert 1997, 1999) show that strong angular 
momentum gradients can result in small, short-lived disk-like 
structures as the flow accretes onto the star. The interaction 
of such an accretion flow with an extant disk warrants further 
investigation. 

Another issue that is unmodeled here is the potential ef- 
fect of protostellar outflows on the accretion flow. There are 
two main effects which could alter the flow; radiation pres- 
sure, and mechanical winds or jets. The influence of stellar 
radiation on BHL accretion has been investigated by Edgar 
& Clarke (2004), who showed that for stellar masses above 
around 10 M Q the accretion rate can be significantly dimin- 
ished, but for lower mass stars like those in this paper radia- 
tive feedback should not present a large impediment to accre- 
tion. The question of jets' and winds' effects on an accretion 
flow was briefly considered in TB08. The ram pressure of a 
jet from a typical T Tauri star could be enough to reverse the 
BHL flow. However, the opening angle of collimated outflows 
tends to be ~ 20° at a distance of 20-50 AU, and smaller at 
larger radii (Hartigan et al. 2004). The solid angle directly 
affected by a jet is then less than about 1 sr at scales com- 
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Figure 5. Snapshots from the simulation with the disk inclined 30° to the inflow, showing the column density on a logarithmic scale. The left column is a 
projection orthogonal to the ISM flow vector, which is from right to left. The right column is the projection along the direction of the ISM flow. Time increases 
downward, as labeled. Only material denser than 2.5 times the background density contributes to the column density. An animation of this simulation can be 
found at http : / /origins . Colorado . edu/~moeckel/BHLpaper /. 



parable to the accretion radius, and assuming that the flow 
intersecting the jet is reversed, this effect is on the order of 
10%. The zero-inclination case we simulated would be more 
greatly affected, as the jet would be directly opposing the ac- 
cretion wake. This would most directly alter the accretion 
rate onto the star, however, and the accretion onto the disk 
would be less impacted. To our knowledge the interaction be- 
tween a collimated outflow and a BHL accretion flow has not 
been modeled numerically, and this presents and interesting 
question for further work, but we estimate that it should be a 
secondary effect. 



With these thoughts about the broader applicability of our 
models in mind, we discuss the main results of our work and 
their implications for disk evolution and planet formation: 

Changes to the disk surface density profile and morphology. 
We have identified two effects that change the surface density 
profile of the disk: accretion onto the disk from the ISM flow, 
and the subsequent redistribution of the extant disk material 
due to absorbing low angular momentum material. While re- 
lated in that the accretion of ISM material leads directly to 
the redistribution of the existing disk, the latter effect is more 
important in terms of changes to the disk surface density pro- 
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Figure 6. Amplitudes of the first eight Fourier components in the azimuth of 
the disk surface density. Top: the face-on case, showing some power in modes 
with m < 5, dominated by m = 1. Bottom: the 30° inclination case, with a 
stronger hi = 1 mode due to the misalignment between the inflow vector and 
disk symmetry axis. Note the difference in scale between the two plots. 



file. With accretion alone the surface density profile would 
increase ~ 10-20% (see figure 3), while with the migration of 
outer material inwards the profile more than doubles across 
much of the disk relative to the initial value. 

A possible consequence of this, pointed out in TB08, is 
the delivery of fresh material to regions of the disk that may 
have been depleted by planet formation, possibly ushering in 
a second era of gas accretion. TB08 note that over several 
Myr accretion can deliver an amount of gas comparable to 
(or even greater than) the original disk mass. This work sug- 
gests that the effect of accretion could be magnified by the 
presence of remnant gas disk, dumping not only the accreted 
ISM gas but the outer disk gas onto the inner disk and any 
planets that may be forming there. The presence of a dissi- 
pative disk can have an important influence on the dynamics 
of multi-planet systems (Moorhead & Adams 2008; Moeckel 
et al. 2008; Thommes et al. 2008; Matsumura et al. 2009), 
and the dynamics of such systems in a disk that is undergoing 
induced global migration may prove interesting. 

A secondary effect is the shrinking outer radius of the disk. 
As long as the disk is accreting ISM material, the inward mi- 
gration of the outer regions continues, as is seen in figure 3. 
Any wispy outer portions of the disk are more susceptible to 
ram pressure, and thus a clean disk edge is maintained. This 
effect would be most pronounced in disks that have undergone 
more accretion, perhaps due to being in a denser-on-average 
region, for instance the center of a cluster potential. As the 
disk redistribution occurs on a timescale comparable to the 
orbital period, this process could take effect quickly and per- 
haps provides a mechanism to generate a relationship between 
disk-radius and the system's radius in a cluster. This would 
work in the same sense as later photoevaporation once a clus- 
ter's O stars begin to dominate their locality, for example as 
in NGC2244 (Balog et al. 2007). 

Triggered gravitational instabilities? With our chosen pa- 
rameters, the disks should be very stable against gravitational 
instabilities that may aid planet formation (e.g. Boss 1997, 
2003; Rice et al. 2006). Despite this, we see spiral features in 
the column density plots (figures 4 and 5). To better under- 
stand these, we examine the azimuthal Fourier amplitudes of 
the surface density, which we compute as 



m-.e 



where mj is the mass of the SPH particle, m is the mode, (fij is 
the azimuthal angle of the particle, and the summation is over 
the particles identified by their angular momentum as part of 
the disk. In figure 6 we plot the amplitude of the first eight 
Fourier modes (normalized by An) for each of the runs, at five 
times each. The Fourier modes for the isolated disk, which 
we do not plot, lack any power above \A m | / \Aq | = 0.003 in the 
first eight modes, confirming that absent the influence of the 
ISM flow, these disks are smooth and stable. In the face-on 
case, the spiral structure begins with a marginal increase in the 
amplitude of the m ~ 4-5 modes which are tracing the spiral 
features most clearly seen in the top panel of figure 4. The low 
amplitudes suggest that these spiral features are not globally 
organized features due to gravitational instability, but rather 
due to a more local reorganization of the disk as the outer 
material is redistributed inwards. 

The more striking feature is the m = 1 mode, which domi- 
nates the Fourier spectrum at later times for = 0° and at all 
times for ; = 30°. This is due to the star accreting a filament of 
material with some momentum perpendicular to the accretion 
axis, so that the star is knocked slightly out of the center of the 
disk. As the disk responds to this off-axis potential, the m = 1 
mode is excited. The signal is much stronger in the inclined 
case, as even accretion of material along the inflow direction 
will move the star slightly out of the disk's axis of symmetry. 
Thus this signal, too, is not due to gravitational instability, and 
the spiral features are unlikely to be signals of disk fragmenta- 
tion leading to planet formation. We note, however, that a disk 
beginning in a more marginal state of stability than ours could 
be pushed into an unstable state by the disk redistribution we 
see. If a factor of ^2 increase in surface density would be 
enough to push a disk into an unstable state, than interactions 
with the ISM like we have modeled here may be sufficient 
to stimulate planet formation by gravitational collapse. For 
the low-mass disk simulated here though, we stress that even 
with the addition of mass from the ISM the Toomre Q param- 
eter remains much larger than unity. The spiral features are 
not due to gravitational instability, but instead the interaction 
of the disk and the star with an accretion flow. This process 
could conceivably masquerade as gravitational instability in 
observations. 

Changes to the disk inclination. Intuitively the addition of 
material with angular momentum misaligned with the rotation 
of the disk seems like a potential route toward the inclination 
difference between the ecliptic plane and the Sun's rotation 
axis. The instantaneous orientation of the disk, measured by 
the average angular momentum of the particles identified as 
belonging to the disk as in the surface density calculation, is 
a bit noisy, but at the end of the simulations when the disk 
has settled down the face-on case exhibits a 0.5° inclination 
change relative to its initial orientation, and the inclined case 
has changed by 1.5°. This mechanism would appear to not 
be enough to account for the Sun's <~ 7° inclination relative 
to the ecliptic plane (Beck & Giles 2005), at least with the 
amount of mass accreted in our simulations. If the accreted 
ISM material overwhelms the initial disk, as shown to be pos- 
sible in TB08, the situation might be different, but a much 
more detailed examination of the timescales for planet for- 
mation versus disk replacement and reorientation would be 
necessary to make this claim. 

To conclude, we have performed simulations studying 
the interaction between a protoplanetary disk and a dense, 
ambient ISM. The accretion rate from the ISM onto the star 
is fairly well described by the standard BHL accretion rate, 
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even in the presence of a disk altering the flow. We have 
showed that these interactions can have a significant effect on 
the physical structure of the disk, with the main effect being 
the redistribution the outer disk inwards. This redistribution, 
and the addition of foreign material to the extant disk, may 
have some implications for planet formation scenarios. 
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